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Introduction 

Two  basic  formulations,  finite-volume  and  finite-difference,  for  the  implemen¬ 
tation  of  high-order  accurate,  essentially  non-oscillatory  (ENO)  shock-capturing 
schemes  have  been  the  subject  of  considerable  interest  in  recent  years.  These 
schemes  achieve  high-order  spatial  accuracy  in  smooth  regions  by  means  of  a  piece- 
wise  polynomial  approximation  operator  that  is  also  designed  to  avoid  oscillations 
associated  with  interpolation  across  steep  gradients.  As  such,  they  are  well  suited 
for  the  study  of  aeroacoustic  and  transition-related  problems  and  may  serve  as  an 
alternative  to  spectral  methods  for  solving  such  problems  when  shocks  or  complex 
geometries  are  involved. 

The  finite-volume  implementation,  first  presented  by  Harten  et  al.1,  is  preferred 
for  its  strict  adherence  to  the  integral  form  in  which  conservation  laws  are  defined. 
The  primary  motivation  for  the  use  of  the  finite-difference  approach  of  Shu  and 
Osher2  is  computational  efficiency.  These  formulations  are  briefly  described,  after 
which  results  of  their  numerical  implementations  are  presented  for  comparison. 
The  intent  of  this  work  is  to  acquaint  the  reader  with  the  relative  merits  of  both 
formulations,  the  circumstances  for  which  each  might  be  useful,  and  some  details  of 
implementation  that  may  be  required  for  a  given  application.  The  performances 
of  both  algorithms  are  compared  for  accuracy,  sensitivity  to  grid  irregularities, 
resolution  of  waves  that  are  oblique  to  the  mesh,  and  computational  efficiency. 

Discrete  Formulations 

The  finite-volume  and  finite-difference  algorithms  compared  in  this  paper  differ 
fundamentally  in  the  way  a  system  of  equations  is  solved.  In  both  cases,  a  weak 
solution  of  a  system  of  conservation  laws  is  ultimately  obtained.  The  conservation 
of  some  quantity  U  in  a  spatial  domain  D  can  be  written 

4  [  UdV  =  -  /  F  n  dS  (1) 

(ft  JD  JdD 

where  F  is  the  flux,  dD  is  the  boundary  of  D,  dV  is  a  volume  element  of  D,  dS  is  an 
element  of  surface  area  on  dD,  and  n  denotes  the  outward  unit  normal  to  dD. 

In  the  finite-volume  approach,  the  conservation  law  itself  in  Eq.  (1)  is  approxi¬ 
mated.  The  spatial  domain  is  discretized,  D  =  { D <},  which  results  in 

= -v,Lfsds  (2) 
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where  Vi  is  the  volume  of  Dt  and 

Ui  =  ^  /  V  dV  (3) 

Vi  JD, 

is  the  cell  average  of  U  in  Dt.  Eq.  (2)  is  solved  for  all  i,  which  yields  a  solution  {(/,} 
of  cell  averages. 

In  the  finite-difference  approach,  a  pointwise  solution  is  desired.  To  this  end,  time 
differentiation  and  spatial  integration  are  interchanged  in  Eq.  (2),  the  divergence 
theorem  is  applied  on  the  right-hand  side  and,  in  the  limit  as  V,  — >  0, 

wu ■  -  -(’■*).  (4) 

Special  care  must  be  taken  when  this  formulation  is  implemented,  because  flux 
conservation  is  not  as  readily  achieved  as  in  Eq.  (2). 

For  integration  in  time,  the  method  of  lines  will  be  employed  for  Eqs.  (2)  and 
(4).  High-order  accurate  Runge-Kutta  methods,  developed  by  Shu  and  Osher2,  are 
implemented  in  the  finite-volume  and  finite-difference  schemes  to  be  compared  in 
this  paper.  Hence,  the  brief  description  of  the  schemes  to  follow  will  concern  only 
the  right-hand  sides  of  Eqs.  (2)  and  (4). 

Both  discrete  algorithms  involve  a  reconstruction  step  followed  by  an  evolution 
step.  What  is  meant  by  reconstruction  is  a  high-order  accurate  polynomial  approx¬ 
imation  at  some  point  in  time.  In  the  finite-volume  formulation,  the  solution  U  is 
reconstructed  from  the  cell  averages  to  high  order  within  each  cell  D,  and  evaluated 
on  the  boundary  dDt.  The  evolution  step  involves  the  solutions  of  the  local  Riemann 
problems  that  arise  from  the  piecewise  continuous  reconstruction  (See,  e.g.,  Refs.  1 
and  3  and  the  references  therein.).  The  spatial  integration  along  the  boundary  dDt  is 
achieved  by  a  correspondingly  high-order  quadrature.  This  method  will  be  referred 
to  throughout  the  paper  as  ENO-FV.  Shu  and  Osher2,4  have  proposed  the  use  of 
the  finite-difference  formulation  in  Eq.  (4),  in  which  the  reconstruction  operator  is 
applied  directly  to  the  pointwise  flux  values.  The  evolution  step  arises  from  a  flux- 
vector  splitting  strategy  that  is  built  into  the  reconstruction.  The  term  ENO-FD 
will  be  used  to  refer  to  this  finite-difference  scheme.  The  pointwise  nature  of  this 
formulation  eliminates  the  need  for  dealing  with  cell  averages  or  the  integration  of 
a  flux  over  the  boundary  of  a  cell.  These  distinctions  between  the  two  formulations 
become  most  important  with  regard  to  the  issue  of  cost,  and  will  be  discussed  in 
more  detail  in  a  later  section. 
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The  most  unique  aspect  of  ENO  reconstruction  operators  is  their  use  of  adaptive 
stenciling.  At  the  i-th  discrete  location,  the  reconstruction  set  that  determines  a 
polynomial  Vi  is  chosen  by  a  searching  algorithm  in  which  decisions  are  based 
on  local  smoothness  criteria.  For  example,  in  one  dimension,  where  the  {D,}  are 
intervals  on  the  Real  line,  Harten  et  al.1  have  suggested  a  hierarchical  algorithm 
that  is  based  on  local  sets  of  A>th  divided  differences  {£*},  where  i  -  m  <  j  <  i  +  m, 
k  =  1,2 , . . . ,  m,  and  m  is  the  desired  degree  of  Vi.  If  a  contiguous  one-dimensional 
stencil  for  Vi  is  defined  by  its  left-most  index  im,  then  the  first-order  choice  is  clearly 
ii  =  i.  Then,  for  k  =  2, 

,t  =  (  ‘*-1  -  1  •  if  K-.I<K‘I  (5) 

^  ijt-l  ,  otherwise 

Because  this  algorithm  allows  the  reconstruction  stencil  to  shift  freely  with  the 
detection  of  any  numerical  gradient,  it  will  be  referred  to  as  “freely  adaptive.” 
Reasons  for  modifying  this  search  will  become  apparent. 

A  few  observations  that  concern  the  available  types  of  reconstruction  operators 
are  in  order.  Within  the  Shu-Osher  approach,  when  conservation  is  desired,  the 
implementation  of  the  high-order  reconstruction  operator  requires  a  uniform  compu¬ 
tational  mesh.  Therefore,  the  application  of  the  ENO-FD  algorithm  on  a  non-uniform 
physical  domain  requires  a  sufficiently  smooth  transformation  to  a  uniform  mesh  if 
third-order  or  higher  accuracy  is  desired.  An  analogous  form  of  this  transformational 
reconstruction  (TR)  can  also  be  implemented  within  the  finite- volume  formulation.  In 
this  case,  the  reconstruction  operator  is  applied  to  a  set  of  volume-weighted  averages 
{ViUi},  and  therefore  sufficient  smoothness  in  a  spatial  transformation  is  required.5 
However,  another  option  exists  for  the  ENO-FV  algorithm.  The  polynomial  approxi¬ 
mation  can  be  performed  in  physical  space,  which  releases  such  burdensome  restric¬ 
tions  on  grid  smoothness.6,7,8  This  latter  procedure  will  be  referred  to  as  physical 
reconstruction  (PR).  Although  this  PR  operator  poses  no  problems  in  one  dimension, 
its  implementation  in  multidimensional  space  can  be  quite  complex  when  the  imple¬ 
mentation  must  allow  for  local  stencil  adaptation.7,8  However,  the  multidimensional, 
finite-volume  TR  operator  can  be  readily  implemented  because  it  is  defined  as  a  prod¬ 
uct  of  one-dimensional  operators.5  When  required,  the  finite-volume  algorithms  will 
be  distinguished  as  ENO-FV-TR  and  ENO-FV-PR. 

One-Dimensional  Rarefaction  Wave 
The  first  test  case  involves  the  solution  of  the  Euler  equations  of  gas  dynamics  in 
one  spatial  dimension,  as  it  pertains  to  the  movement  of  a  right-traveling  rarefaction 
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wave  through  a  domain  of  highly  varying  mesh  spacing.  The  effects  of  two  mesh 
transformations  will  be  examined.  One  transformation  has  the  required  smoothness 
for  a  fourth-order  TR  operator  and  the  other  does  not. 

For  both  transformations,  the  uniform  computational  domain  is  given  by 
{ -6  <  £  <  6 } .  The  subset  { -4  <  £  <  4 }  is  divided  into  five  intervals: 
[-4,  -3],  [-3,  -1],  [-1, 1],  [1,3],  and  [3,4] .  In  each  of  these  interior  intervals,  the  map¬ 
ping  x  =  x(£)  causes  the  physical  mesh  width  Ax  to  vary  rapidly;  the  intervals 
{  -6  <  £  <  -4 }  and  { 4  <  £  <  6 }  are  mapped  uniformly. 

In  the  A:-th  interior  interval,  the  first  transformation  is  of  the  form 

*({)  =  &  +  sin  ^  (£-&■)  (6) 

where  £*  is  an  element  of  the  set  { -4,  -2, 0, 2, 4} .  The  mesh  spacing  Ax  in  the  uniform 
regions  (fore  and  aft)  is  determined  such  that  the  connections  at  £  =  ±4  are  smooth. 

To  generate  a  smoother  grid  on  {-4<£<4},  a  mapping  of  the  following  form 
is  used: 

*(£)  =  £+  ^  y  sin  (tt£)  -  ^  sin  (2tt£)  +  ^-  sin  (3?r£)  (7) 

The  parameters  a  and  0  are  determined  so  that  the  ratio  of  the  maximum  to 
minimum  values  of  are  identical  for  the  two  grids,  and  that  the  physical  distance 
between  x(-4)  and  x(4)  is  the  same  for  both. 

Fig.  1(a)  illustrates  the  similar  behavior  of  these  two  transformations  on 
0  <  £  <  2 .  The  seemingly  odd  formula  in  Eq.  (7)  was  chosen  because  its  derivative 
is  of  the  form 

X{=  a-  f3  sin6  (|  £  ) 

The  transformation  derivatives  are  plotted  on  the  same  interval  in  Fig.  1(b).  The 
value  x$  represents  the  derivative  normalized  by  its  maximum  value.  For  the 
transformation  (6),  x^  is  discontinuous  for  £  =  ±1,±3,  which  makes  for  an  overall 
C1  grid.  The  mesh  produced  by  Eq.  (7)  will  be  C°°  on  the  nonuniform  region  and  C6 
overall,  because  of  the  connections  with  the  uniform  intervals  at  £  =  ±4  . 

The  initial  solution  consists  of  an  isentropic  expansion,  smoothly  distributed  on 
-6.5  <  x  <  —5.5,  with  constant  states,  Kj  and  Ui,  on  the  right  and  left,  respectively. 
The  strength  of  this  rarefaction  is  determined  by  requiring  a  mean  temperature 
change  of  ±5  percent  across  the  wave  and  supersonic  flow  on  either  side.  The  problem 
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is  nondimensionalized  with  respect  to  the  mean  solution.  For  t  >  0 ,  the  rarefaction 
wave  moves  to  the  right,  and  its  behavior  is  monitored  within  the  nonuniform  region 
of  the  mesh  until  t  =  4.0. 

Both  algorithms  used  here  are  fourth-order  accurate  in  the  L\  sense.  Fig.  2 
compares  solutions  for  the  two  formulations  by  x-t  contour  plots  of  the  density,  for 
0.0  <t  <  4.0 ,  on  a  mesh  of  120  intervals  with  the  C1  transformation  given  by  Eq.  (6). 
The  freely  adaptive  stencil  algorithm  in  Eq.  (5)  is  employed.  All  of  the  contour  plots 
in  this  section  represent  the  entire  solution  range,  p\  <  p  <  p2,  with  20  equally  spaced 
levels.  The  ENO-FD  solution  exhibits  significant  qualitative  error  in  neighborhoods 
about  the  second-derivative  discontinuities  in  the  mesh.  In  Fig.  2(b),  such  problems 
are  not  apparent  when  the  solution  is  computed  with  the  freely  adaptive  ENO-FV- 
PR  algorithm. 

Fig.  3  depicts  the  use  of  the  same  fourth-order  ENO-FD  algorithm  on  the  C6 
grid  generated  by  Eq.  (7).  With  the  required  mesh  smoothness  to  support  the 
reconstruction,  no  visible  distortions  appear  as  on  the  C1  grid.  However,  the  fourth- 
order  design  accuracy  is  not  achieved  for  this  solution  when  density  error  is  measured 
with  respect  to  the  Lj  norm  in  a  region  to  the  left  of  the  rarefaction  at  Z  =  4.0.  The 
same  accuracy  problem  was  found  with  the  freely  adaptive  ENO-FV-PR  scheme.  The 
results  of  both  mesh-refinement  studies  are  shown  on  a  log-log  plot  in  Fig.  4(a).  The 
spatial  discretizations  employed  in  this  study  are  120,  240,  480,  and  960  intervals. 
The  numbers  on  the  loci  represent  the  computational  order  of  accuracy  as  measured 
between  the  two  finest  meshes. 

A  similar  loss-of-accuracy  phenomenon  was  reported  by  Rogerson  and  Meiberg9, 
which  prompted  a  response  from  Shu10  that  the  problem  arises  by  allowing  the 
stencil  to  adapt  too  freely.  Shu  has  suggested  that  the  stencil  adaptation  algorithm 
be  modified  to  bias  the  stencil  towards  one  that  is  stable,  in  the  sense  of  linear 
stability  analysis.  In  the  present  application,  the  resulting  stencil  is  one  which  is 
upwind  biased.  This  biasing  can  be  done  by  implementing  a  factor  a  in  the  stencil 
search  in  Eq.  (5),  vis. 

i*  *  (**->~ 1  ’  if  <^-il<*R|*iLl  (9) 

l,  ik-i  ,  otherwise 

where  (<7£,  ctr)  =  (1,<t)  or  (a,  1) ,  for  biasing  to  the  left  or  right,  respectively,  with 
a  >  1 . 

Fig.  4(b)  shows  grid-refinement  results  that  are  analogous  to  those  in  Fig.  4(a) 
using  the  modification  in  Eq.  (8)  with  a  =  2.0;  however,  a  fourth-order  error  reduction 


is  still  not  evident  because  the  error  is  measured  in  a  region  where  numerical 
gradients  are  extremely  small  and  ratios  of  neighboring  gradients  may  be  much 
larger  than  the  chosen  biasing  parameter.  Atkins11  has  suggested  that  an  additional 
constraint  on  the  stencil  adaptation  is  required.  The  purpose  of  this  constraint  is 
to  bias  the  stencil  toward  one  that  is  stable  wherever  the  solution  is  smooth.  A 
parameter  installed  for  this  purpose  can  be  considered  a  lower  threshold  for  the 
magnitudes  of  the  local  differences  {&*■} ,  below  which  a  stencil  is  forced  toward 
a  stable  target,  regardless  of  the  relative  magnitudes  of  neighboring  numerical 
gradients.  The  present  implementation  of  such  a  parameter  can  be  written 

if  &ik-i  <  t  and  <  e 

then  ik=  H  ■>  (10) 

where  e  is  a  small  parameter,  and  i*k  defines  the  stable  stencil  at  the  k- th  level. 

Another  grid-refinement  study  was  performed,  with  Eqs.  (8)  and  (9)  wherein 
a  =  2.0  and  e  =  0.01.  The  results  are  shown  in  Fig.  5(a)  and  are  much  more  consistent 
with  the  scheme’s  design  accuracy.  In  fact,  the  ENO-FV-PR  algorithm,  coupled  with 
the  stenciling  modifications  in  Eqs.  (8)  and  (9),  yields  fourth-order  computational  ac¬ 
curacy  on  the  C1  grid,  which  is  shown  in  Fig.  5(b).  The  second-order  convergence 
exhibited  by  the  ENO-FD  algorithm  in  this  plot  is  expected.  However,  Fig.  6  shows 
that  the  qualitative  error  in  the  ENO-FD  solution  on  the  C1  grid,  that  was  previ¬ 
ously  shown  in  Fig.  2(a),  has  been  significantly  reduced  with  these  stencil  biasing 
modifications. 

At  this  point,  the  apparently  greater  robustness  of  the  ENO-FV  algorithm  can  be 
attributed  to  the  fact  that  this  formulation  allows  for  a  physical  reconstruction  that 
is  not  available  for  the  finite-difference  algorithm.  However,  as  previously  noted, 
multidimensional  extensions  for  this  more  generalized  operator  can  be  complicated 
and  costly.  Therefore,  some  results  produced  with  the  ENO-FV-TR  algorithm  will  be 
discussed.  In  Fig.  7(a),  when  the  freely  adaptive  version  of  this  algorithm  is  imple¬ 
mented  on  the  C6  grid,  significant  error  exists  that  was  not  evident  in  the  analogous 
case  for  the  ENO-FD  solution  in  Fig.  3.  This  error  is  due  to  the  reconstruction  of 
volume-weighted  averages,  in  which  a  rapidly  varying  mesh  will  have  an  inordinate 
effect  on  the  stencil  choice  in  a  region  of  small  solution  gradients.  However,  biasing 
the  reconstruction  stencils  with  Eqs.  (8)  and  (9)  is  not  enough  to  entirely  rid  the 
solution  of  these  grid  distortion  errors,  as  shown  in  Fig.  7(b).  The  problem  is  still 
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related  to  the  reconstruction  of  volume-weighted  averages,  in  particular  to  the  fact 
that  a  uniform  flow  is  not  preserved  by  this  procedure  if  the  analytic  transformation 
is  used  for  the  necessary  discrete  values  {x^}  .  Therefore,  instead  of  using  Eq.  (7b)  to 
compute  the  derivatives,  these  derivatives  are  numerically  approximated  in  a  man¬ 
ner  that  enables  the  reconstruction  operator  to  preserve  the  free  stream  exactly5. 
The  resulting  solution,  with  this  final  modification,  is  shown  in  Fig.  7(c).  In  this 
final  form,  the  numerical  accuracy  of  the  ENO-FV-TR  algorithm  is  found  to  perform 
to  design  as  shown  in  Fig.  8(a).  The  ENO-FD  and  ENO-FV-PR  results  from  Fig.  5(a) 
are  repeated  in  Fig.  8(b),  for  comparison. 

Two-Dimensional  Channel  Flow 

The  ENO-FD  and  ENO-FV  algorithms  are  now  compared  in  two  spatial  dimen¬ 
sions.  The  test  case  involves  a  steady,  subsonic  flow  in  a  channel  of  varying  area. 
Although  high-order  ENO  schemes  are  clearly  designed  with  unsteady  solutions  in 
mind,  many  such  solutions  of  interest  can  be  considered  as  the  imposition  of  per¬ 
turbations  upon  a  steady  flow.  It  is  therefore  essential  that  a  steady  flow  be  accu¬ 
rately  predicted  in  order  to  obtain  meaningful  results  from  unsteady  problems  of  an 
aeroacoustic  or  transitional  nature.  Unless  otherwise  stated,  the  remainder  of  the 
applications  of  the  ENO-FV  algorithm  will  involve  the  TR  type  of  reconstruction. 

This  channel  flow  solution  is  assumed  to  be  governed  by  the  two-dimensional 
Euler  equations  and  is  computed  with  both  the  ENO-FD  and  ENO-FV  algorithms 
on  two  different  geometries.  For  the  two  channels  under  consideration,  the  length- 
to-height  ratio  is  L/H  =  1.5,  with  constant-area  sections  fore  and  aft  of  a  section 
of  varying  area.  Each  constant-area  section  has  length  L/ 5;  at  the  throat,  the  con¬ 
striction  is  10  percent  of  H.  The  most  significant  difference  between  the  two  geome¬ 
tries  is  the  order  of  smoothness  to  which  the  constant-area  sections  are  connected 
to  the  center  section.  Both  channels  can  be  described  as  follows.  Let  the  rectangle 
{0  <  £  <  L)  x  {0  <  r}  <  H}  denote  the  computational  domain  for  each  channel.  The 
identity  maps  (£,  77)  to  the  physical  point  (x,  y) ,  for  £  in  [0,  L/5]  U  [4L/5,  L]  and  all  77. 
On  the  interval  {1/5  <  (  <  4L/5} ,  the  geometry  for  the  varying-area  section  is  given 
by  a  transformation  of  the  form 


x  =  (11a) 

#  =  £««)  +  (  1 -•£)»!«)  (lib) 

where  2/1  (0  and  y2(0  are,  respectively,  the  equations  for  the  bottom  and  top  walls. 
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The  first  middle  section  considered  is  determined  by  walls  of  the  form 

^•(0  =  sin4  (6,£  +  C{)  +  di  ,  *  =  1,2  (12) 

where  the  values  of  the  coefficients  are  set  for  the  desired  throat  constriction  and 
a  smooth  connection  to  the  constant-area  sections.  Three  continuous  £  derivatives 
of  y  exist  at  £  =  I/5,4L/5 ,  for  all  tj.  Therefore,  Eqs.  (10)  and  (11)  generate  a  C3 
geometry.  Fig.  9  illustrates  this  geometry  on  a  60x40  mesh.  The  second  channel 
geometry  differs  from  the  first  only  in  that  the  walls  in  the  middle  section  are  given 
by  polynomials  of  the  form 


yi(t)  =  ai  +  bit2  +  ci?  ,  *  =  1,2  (13) 

In  this  case,  the  connections  of  the  middle  to  the  outer  sections  are  continuous  in  £ 
to  only  one  derivative,  for  0  <  tj  <  H . 

The  desired  solution  is  that  of  a  steady-state  flow  that  is  caused  by  a  uniform, 
parallel  free  stream  entering  the  channel  at  x  =  0 .  This  solution  is  achieved  by 
solving  the  time-dependent  Euler  equations  with  the  fourth-order  algorithms.  At 
t  =  0 ,  the  solution  in  the  entire  channel  is  set  to  free  stream  with  a  Mach  number  of 
0.3.  Tangency  is  imposed  on  the  walls,  and  a  non-reflecting  boundary  condition12  is 
applied  at  the  inflow  and  outflow.  The  stencil-biasing  modifications  in  Eqs.  (8)  and  (9) 
are  essential  for  convergence  of  the  solution  to  a  steady  state,  and  are  implemented 
with  a  =  2.0  and  e  =  0.01 .  Flux  residuals  were  readily  driven  to  machine  zero  in 
both  test  cases. 

These  solutions  were  computed  by  both  algorithms  on  a  sequence  of  successively 
refined  grids,  and  the  solutions’  errors  were  determined  by  deviation  from  isentropy, 
as  measured  by  the  quantity  S  =  P/p'1 .  The  refinement  sequence  employed  here 
is  15x10,  30x20,  45x30,  60x40,  and  90x60.  The  results  of  the  grid-refinement 
study  on  the  C3  geometry  are  shown  in  Figs.  10(a),  (b),and  (c).  The  “Global”  er¬ 
ror  is  computed  over  the  entire  computational  domain,  the  “Wall”  error  is  com¬ 
puted  only  at  the  points  along  one  wall,  and  the  “Interior”  error  is  computed  on 
{L/ 4  <  (  <  3L/4}  U  {H/ 4  <  17  <  3///4} .  Both  algorithms  evidently  perform  at  or  near 
design  accuracy,  as  expected  for  fourth-order  schemes  on  a  C3  mesh. 

In  Figs.  10(d),  (e),  and  (0  are  the  grid  refinement  results  for  the  C1  mesh. 
As  expected,  second-order  results  are  obtained  on  the  wall  for  both  algorithms. 
However,  the  finite-volume  algorithm  performs  at  third-order  accuracy  with  respect 
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to  the  global  error  and  at  design  accuracy  on  the  interior,  while  the  finite-difference 
algorithm  shows  second-order  accuracy  in  all  three  measures.  These  results  suggest 
that,  for  the  ENO-FD  algorithm,  the  second-order  entropy  error  that  arises  from 
the  non-smoothness  at  the  section  connections  is  propagating  into  the  interior.  This 
assumption  is  supported  by  Fig.  11,  in  which  the  quantity  log  S  is  plotted  along  the 
center  line  of  a  60x40  mesh.  Clearly,  in  the  finite-difference  solution,  there  is  a 
lower-order  entropy  error  within  the  middle  section  than  exists  in  the  constant-area 
sections. 

Oblique  Sod’s  Problem 

The  final  test  under  consideration  is  selected  to  compare  the  capabilities  of  the 
two  algorithms  to  resolve  waves  that  are  oblique  to  the  computational  mesh.  Sod’s 
problem13  will  be  solved  in  two-dimensional  space  so  that  the  planar  waves  produced 
will  propagate  at  various  angles  of  incidence  with  respect  to  a  rectangular  grid.  The 
intent  is  not  only  to  inspect  the  qualitative  resolution  of  the  oblique  waves,  but  also 
to  quantify  the  manner  in  which  each  algorithm  detects  an  oblique  wave  with  respect 
to  its  detection  of  a  wave  that  is  normal  to  the  me^n. 

Sod’s  problem  is  a  Riemann  problem  which  is  subject  to  the  solution  of  the  Euler 
equations  and  determined  by  initial  conditions  that  consist  of  a  thermodynamic 
discontinuity  imposed  upon  a  fluid  at  rest.  The  magnitudes  of  the  density  and 
pressure  jumps  are,  from  left  to  right,  pl/pr  =  8  and  Pl/Pr  =  10,  respectively. 
Such  an  initial  solution  is  installed  on  a  two-dimensional  Cartesian  grid  so  that  the 
discontinuity  is  a  straight  line  that  makes  an  angle  0  with  the  x  axis.  The  physical 
mesh  varies  with  0  in  the  following  manner.  A  computational  mesh  is  defined  on 
the  rectangle  [0,1]  x  [0,  H] ,  with  a  length-to-height  ratio  L/H  =  6.  If  0  <  0  <  tt/2, 
then  the  physical  domain  is  {0  <  x  <  Lg}  x  {0  <  y  <  Hg} ,  and  is  determined  by 
Hg  =  H/sinO  and  Lg  =  LI  sin  0 .  This  scaling  achieves  the  same  grid  resolution 
normal  to  the  wave  propagation  on  a  given  mesh  at  some  fixed  time  for  all  choices 
of  0.  The  mesh  discretization  in  each  direction  is  chosen  such  that  Ax  =  A y . 

At  t  =  0,  the  initial  discontinuity  is  positioned  at  (x,y)  -  (31/8,0)  and  inclined 
at  the  angle  0.  The  angles  of  inclination  are  chosen  so  that  tan  0  is  an  integer  and 
that  the  boundary  conditions  at  y  =  0  and  y  =  Hg  can  be  determined  in  a  “shifted- 
periodic”  manner.  In  particular,  the  test  angles  are  0  =  arctan  1 ,  arctan  2,  and  arctan  4 , 
in  addition  to  the  one-dimensional  problem  0  =  tt/2  .  The  solution  is  computed  to 
t  =  1.2  The  stencil  biasing  in  Eq.  (8)  is  used  in  the  reconstruction  procedures  with 
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a  =  2.0. 


Fig.  12(a)  represents  the  solution  at  t  =  1.2,  on  a  96x16  grid,  with  0  =  arctan  1, 
with  the  use  of  the  ENO-FV  algorithm.  The  three  wave  structures  are,  from  left  to 
right,  an  expansion  wave,  a  contact  discontinuity,  and  a  shock.  The  axis  variables  x' 
and  y'  represent  the  physical  coordinates  scaled  by  sin  0 .  Fig.  12(b)  depicts  the  same 
solution  as  in  Fig.  12(a),  for  both  algorithms,  along  the  grid  line  y'  =  0 .  On  this  level, 
the  qualitative  difference  between  the  solutions  is  barely  detectable.  Therefore,  in 
Fig.  13,  the  solutions  are  examined  in  the  vicinity  of  the  contact  discontinuity  for 
all  test  values  of  8.  (The  successive  loci  in  this  plot  have  been  shifted  upward  and 
to  the  right.)  Although  the  ENO-FV  scheme  produces  a  slightly  steeper  numerical 
gradient  at  the  location  of  the  discontinuity,  such  a  small  difference  can  be  judged 
insignificant.  Fig.  13,  then,  suggests  that  both  algorithms  perform  equally  well 
with  respect  to  oblique  waves  irrespective  of  6.  On  a  more  quantitative  level, 
the  differences  between  the  oblique  cases  and  the  one-dimensional  case,  for  each 
algorithm,  have  been  measured.  These  results  are  plotted  in  Figs.  14(a)  and  14(b). 
By  this  measure,  the  finite-volume  scheme  performs  marginally  better,  particularly 
with  respect  to  the  rarefaction  wave. 

Cost  Comparison 

As  with  any  numerical  algorithm,  the  cost  of  implementing  either  the  ENO-FD 
or  ENO-FV  scheme  is  a  major  concern.  Conceptually,  the  two  formulations  can  be 
made  equally  cost  effective  for  high-order  solutions  of  one  dimensional  problems  or  for 
first  or  second-order  accurate  solutions  in  multiple  dimensions.  However,  for  third 
or  higher-order  accuracy  in  two  or  more  dimensions,  the  algorithms  are  radically 
different.  This  difference  translates  to  a  significant  disparity  in  cost. 

With  regard  to  cost,  the  ENO-FD  algorithm  has  a  clear  advantage  when  applied 
to  multidimensional  problems.  This  advantage  can  be  entirely  attributed  to  the  fact 
that  the  finite-difference  operator  solves  a  system  of  equations  in  a  pointwise  man¬ 
ner.  In  this  case,  a  high-order  multidimensional  reconstruction  can  be  accomplished 
in  a  dimension-by-dimension  fashion.  If  the  algorithm  is  r-th  order  accurate  in  the 
L\  sense,  then  the  reconstruction  stencil  in  k  dimensions  contains  k(r  -  1)  -f  1  points. 
The  dimension-by-dimension  approach  to  reconstruction  is  not  possible  within  the 
finite-volume  formulation,  if  third-order  or  better  accuracy  is  required.  Because  the 
solution  at  any  time  is  in  the  form  of  cell  averages,  the  multidimensional  ENO-FV- 
TR  reconstruction  is  implemented  as  a  product  of  one-dimensional  operators,  and 


10 


thereby  requires  r*  cells  in  a  reconstruction  stencil.  Moreover,  the  Gaussian  quadra¬ 
ture  used  to  integrate  the  flux  on  each  cell  boundary  requires  multiple  solutions 
of  the  Riemann  problem  on  each  cell  interface.  The  number  of  points  nq  required 
in  this  quadrature  on  an  interface  is  nq  =  [int(  (r+l)/2)]*-1  ,  where  int(-)  denotes 
integer  truncation.  This  high-order  quadrature  is  entirely  avoided  within  the  finite- 
difference  approach.  It  should  be  noted  that  both  formulations  require  the  same 
amount  of  logic  in  the  stencil  adaptation,  which  is  the  most  expensive  part  of  either 
algorithm. 

Fig.  15  depicts  the  relative  CPU  cost  of  the  two  algorithms,  for  solutions  of 
the  Exiler  equations  in  two  dimensions,  for  2nd,  3rd,  4th,  and  5th-order  accurate 
applications.  CPU  times  were  measured  on  a  Cray-YMP  and  normalized  with 
respect  to  the  2nd-order  ENO-FD  cost.  An  estimated  cost  comparison  for  three- 
dimensional  problems  can  be  obtained  by  the  following  reasoning.  For  a  specified 
order  of  accuracy,  let  C\  denote  the  cost  of  solving  one  equation  at  a  point  in  one 
dimension;  this  is  the  same  for  either  algorithm.  Then,  because  of  its  dimension-by- 
dimension  implementation,  the  estimated  cost  incurred  by  the  ENO-FD  algorithm 
in  solving  the  Eider  equations  in  k  dimensions  on  a  grid  of  Nk  points  is 

CFD  =  Ci{k  +  2)kNk  (14) 

This  formula  predicts  a  cost-per-point  increase  factor  of  8/3  when  the  one-dimensional 
ENO-FD  algorithm  is  extended  to  two  dimensions.  This  value  has  been  supported 
by  computer  measurements.  The  corresponding  factor  for  the  extension  from  two  to 
three  dimensions  is  15/8.  The  estimation  of  a  three-dimensional  cost  for  the  ENO- 
FV  algorithm  cannot  be  done  by  such  a  simple  linear  extrapolation.  There  are  two 
significant  costs  in  addition  to  the  base  value  given  by  Eq.  (13).  One  is  the  extension 
of  the  finite-volume  reconstruction  to  a  three-dimensional  product.  The  other  arises 
from  the  additional  3  (n9  -  1)  flux  computations  that  are  required  in  the  integration 
over  the  surface  of  each  cell.  These  additional  costs  were  extrapolated  from  two- 
dimensional  CPU  measurements  on  a  Cray-YMP.  Fig.  16  represents  the  estimated 
three-dimensional  cost  comparison  where,  again,  the  data  are  normalized  by  the 
2nd-order  cost. 

Concluding  Remarks 

For  accuracy  on  a  sufficiently  smooth  mesh  and  resolution  of  oblique  waves,  both 
algorithms  perform  equally  well.  The  finite-volume  implementation  is  less  sensitive 
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to  derivative  discontinuities,  whether  in  the  computational  mesh  or  in  the  solution. 
In  particular,  the  ENO-FV-PR  algorithm  has  the  capacity  to  perform  at  design 
accuracy,  independent  of  the  mesh.  Although  the  generalized  multidimensional 
adaptive-stenciling  implementation  of  this  algorithm  has  not  sufficiently  matured, 
some  promising  results  can  be  found  in  Ref.  8.  Either  of  the  multidimensional  finite- 
volume  algorithms  is  significantly  more  costly  than  the  finite-difference  algorithm. 
Therefore,  if,  for  a  given  application,  the  computational  domain  is  known  to  be 
sufficiently  smooth  and  can  be  suitably  structured,  then  the  ENO-FD  algorithm  is 
the  method  of  choice.  However,  for  problems  with  complex  geometries,  it  might  pay 
to  use  the  more  expensive  algorithm  if  the  grid  is  significantly  less  costly  to  generate 
in  a  less  restrictive  fashion. 
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Fig.  4.  Lj  density  error,  t  =  4.0,  C4  grid,  (a)  Freely  adaptive  stencil,  (b)  Biased  stencil  ( <t  -  2.0) 
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Fig.  8.  L\  density  error,  t  =  4.0,  C6  grid,  (a)  Biased  stencil,  free-stream  preserving,  (b)  Biased  stencil 

Or  =  2.0,  e  =  0.01). 
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Fig.  9.  The  C4  geometry  of  Eqs.  (10)  ami  (11). 
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fig.  It.  Li  entropy  error*,  steady  2D  channel  flow. 
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Fig.  12.  (a)  Density  contours,  ENO-FV,  9  =  arctan  1.  (b)  Density,  y'  =  0.0,  0  =  arctan  1 
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Fig.  13.  Density,  contact  discontinuity,  y'  =  0.0. 
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Fig.  15.  Relative  CPU  cost  of  the  two  algorithms  in  two  spatial  dimensions. 


Fig.  16.  Relative  CPU  cost  of  the  two  algorithms  in  three  spatial  dimensions  (estimated). 
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